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Abstract. Rhythmic behaviors in neural systems often combine features of limit cycle dynam- 
ics (stability and periodicity) with features of near heteroclinic or near homoclinic cycle dynamics 
(extended dwell times in localized regions of phase space) . Proximity of a limit cycle to one or more 
saddle equilibria can have a profound effect on the timing of trajectory components and response to 
both fast and slow perturbations, providing a possible mechanism for adaptive control of rhythmic 
motions. Reyn showed that for a planar dynamical system with a stable heteroclinic cycle (or sep- 
aratrbc polygon), small perturbations satisfying a net inflow condition will generically give rise to a 
stable limit cycle (Rcyn, 1980; Guckenheimer and Holmes, 1983). Here we consider the asymptotic 
behavior of the infinitesimal phase response curve (iPRC) for examples of two systems satisfying 
Reyn's inflow criterion, (i) a smooth system with a chain of four iiyperbolic saddle points and (ii) 
a piecewise linear system corresponding to local linearization of the smooth system about its saddle 
points. For system (ii), we obtain exact expressions for the limit cycle and the iPRC as a function 
of a parameter fi > representing the distance from a heteroclinic bifurcation point. In the ^ — > 
limit, we find that perturbations parallel to the unstable eigenvector direction in a piecewise linear 
region lead to divergent phase response, as previously observed (Brown, Moehlis and Holmes (2004) 
Neural Computation) . In contrast to previous work, we find that perturbations parallel to the stable 
eigenvector direction can lead to either divergent or convergent phase response, depending on the 
phase at which the perturbation occurs. In the smooth system (i), we show numerical evidence of 
qualitatively similar phase specific sensitivity to perturbation. 
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1. Introduction. Animals often generate specific sequences of behavior, for ex- 
ample the movements of the limbs during walking, the feeding apparatus while chew- 
ing and swallowing, or body undulations in swimming. When a repeated sequence of 
motions can be produced reliably, the pattern generator circuit controlling the behav- 
ior is typically modeled as an autonomous system of ordinary differential equations 
admitting a stable isolated periodic orbit, i.e. a limit-cycle oscillator [30, 57]. Limit 
cycle oscillators have played a fundamental role in understanding the generation and 
control of repetitive motions underlying swimming in the lamprey [12, 15], neural ac- 
tivity and bursting [17, 31], the gaits of quadrupeds [13, 27] bipeds [35] and monopeds 
[9, 14], as well as cardiac and respiratory activity [19, 42, 43], inter alia. 

Another class of systems generating reproducible sequences of activity has been 



proposed under the rubric of stable heteroclinic sequences [5] or stable heteroclinic 
channels [37]. A dynamical system possesses a heteroclinic sequence if there exists 
a chain of hyperbolic saddle fixed points for which the imstable manifold of each 
saddle intersects the stable manifold of the next [2]; they generalize the notion of a 
stable heteroclinic cycle, an attractor comprising a finite collection of saddle points 
with heteroclinic connections linking them in a repeating chain [7, 28, 48, 49]. A 
heteroclinic cycle is stable or attracting when the product around the cycle of the 
saddle values is strictly greater than unity; the saddle value for a hyperbolic saddle 
point with eigenvalues A„ > > > Ag > • • • > A^~^ is the ratio — A^/A„. Stable 
heteroclinic sequences and cycles have been proposed to provide a framework within 
nonlinear dynamics for understanding a range of phenomena, including olfactory pro- 
cessing in insects [36], search behavior in the marine mollusk Clione limacina [56], 
"winnerless competition" in neural circuits [1, 4] and ecological models [2], genesis of 
network-dependent bursting activity [34], and the balance of emotion and cognition 
in behavioral control [3] as well as other areas [38]. 

Because they require the intersection of one dimensional unstable and codimen- 
sion one stable manifolds, the saddle connections comprising a heteroclinic cycle are 
structurally unstable [6]. For planar systems, Reyn showed that a phase portrait 
containing a separatrix polygon generically forms a limit cycle when subject to a per- 
turbation satisfying a net inflow condition [41], provided the unperturbed heteroclinic 
cycle is attracting. Trajectories near an unperturbed attracting heteroclinic cycle 
traverse the cycle with longer and longer return times; the trajectory along the cycle 
itself has "infinite period" . 

As an example, consider the system defined on the 2-torus < yi < 27r, i £ {1, 2}. 

f{yi,y2) = cos(2/i) sin(2/2) -Fasin(22/i) (1.1a) 
g{yi,y2) = -sin(yi)cos(y2) +Q!sin(2y2) (1-lb) 

As shown in Figure 1.1a, this system has four saddle points connected by heteroclinic 
connections. The saddles have identical eigenvalues A„ = 1 — 2a, Ag = — 1 — 2q! for 

a saddle value oi V = (j^) . If a e (0,1/2) then V > 1, indicating that the 
heteroclinic 4-cycle is attracting. 

Although the heteroclinic cycle connecting the four saddle points is structurally 
unstable, any perturbation of the vector field that pushes the unstable manifold of 
each saddle towards the interior of the cycle relative to the stable manifold of the 
next saddle will generically lead to the formation of a stable limit cycle. For example, 
consider the effects of "rotation" of the flow of Equations 1.1 parametrized by 1 > 
At > 0. 

^ = f{yi,y2) + fj'9{yi,y2) (i.2a) 

^ = 5(2/1,2/2) - m/(2/i, 2/2) (l-2b) 

We will refer to this system as the sm,ooth system,, in contrast to the piecewise linear 
system to be introduced subsequently. As shown in Figure 1.1b, for n > the 
heteroclinic connections are broken, but a limit cycle has been created inside of the 
saddles that passes near each of the saddle points. 
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Fig. 1.1: Parametric perturbation of an attracting heteroclinic cycle to a stable limit 
cycle. Note solid line is x-nuUcline and dashed is y-nuUcline. (a) A trajectory of the 
smooth toroidal system given by Equations 1.1, passing near four distinct saddles each 
with eigenvalues = — 1 — 2a and Au = 1 — 2a. The trajectory shown begins near the 
unstable spiral point at location (tt, tt) near the center of the plot. The trajectory was 
integrated for a total time of 200 units using a = 7/30 and fj, = . It passes closer to 
each successive saddle, slowing progressively, (b) The perturbed system, Equations 
1.2, when a = 7/30 and /i = 0.01, forms a stable limit cycle passing close to the 
four saddle points. Please see corresponding movie (file smooth. mpg); this animation 
cycles through phase portraits for the smooth system for values of fi varying from 
to 1 /2 and back again. At /i = 2a the limit cycle collapses via a Hopf bifurcation to 
a single stable fixed point located between the four surrounding saddles. 



For any small positive value of /i, the system 1.2 will have a stable limit cycle, 
7(i) — ■y{t + T), with period T{fi) — >• oo as /x — >■ 0+. For ^ > we may identify a 
phase 9 G [0, 0max.) with each point on the limit cycle by chosing an arbitrary point 
7o to have phase ^(70) = and requiring d9/dt = 2n /T{fi). Typically one chooses 
^max to equal either 2tt or unity. Because of the fourfold symmetry of the systems 
considered here it will be convenient to set ^max = 4 throughout, so that traversal of 
each quarter of the limit cycle will correspond to a unit increment in phase. Figure 
1.2 illustrates the time course of trajectories of Equations 1.2 for a = 7/30 and 
M e {10-3,0.1,0.3,0.45}. 

To each point ^0 = {^OiVq) in the basin of attraction of the limit cycle we may 
assign an asymptotic phase 4'{£,o) G [0,4) satisfying — 7(t — 9{£_o)T/'i)\\ — )• as 
t — >■ 00, where £^{t) and 7(i) are solutions of 1.2 satisfying initial conditions ^(0) = ^0 
and 7(0) — 7o respectively, and ||(a;,y)|| = |a:;| + \y\. The level curves of </>, called 
isochrons, foliate the basin of attraction; see [31] for examples. 

An important question for understanding the dynamics of control in central pat- 
tern generators is how such systems combine robustness to perturbation with flexibil- 
ity to adapt behavior to variable environmental or physiological conditions.^ Many 

^Whether originating endogenously (neural noise, internal control mechanisms), or exogenously 
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Fig. 1.2: Time plots of limit cycle trajectories of the smooth system with various 
values of /i. When /i <C 1 the time plots exhibit prolonged dwell times due to slow 
transits past the saddle points and relatively rapid transits between them. As /U in- 
creases, the limit cycle moves away from the saddle points and the speed becomes 
more uniform (note the change in horizontal scale). As a result, the time plots be- 
come more sinusoidal resembling an Andronov-Hopf oscillator. The first and second 
ordinates are shown in black and gray, respectively, with a. = Prom top to 

bottom n = 10-^,0.1,0.3,0.45. Compare Figure 2.3. 



biological systems combining repetitive behavior and behavioral or metabolic control 
exhibit limit cycle behavior that is strongly influenced by passage of trajectories near 



(environmental fluctuations), perturbations of the system's dynamics occur on a variety of time scales. 
For clarity we will limit discussion to two extremes of fast and slow perturbations. Perturbations 
occurring on slow time scales relative to other system dynamics will be referred to below as static 
or parametric perturbations, for example fixing a value /i > in Equations 1.2. Fast perturbations 
will be approximated as instantaneous trajectory dislocations. 
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one or more unstable fixed points or quasioquilibria. A recent model for the genera- 
tion of multiple rhythmic states in a respiratory CPG featured slow excursions along 
equilibrium surfaces of model neurons voltage dynamics interspersed with fast jumps 
between these surfaces [42] . In this model, control of the rhythm through changes in 
the period as well as the duration of different functional phases (inspiratory phase, 
expiratory phase) could be effected by modulatory signals making small changes to 
the dynamics of escape and release from inhibition, thereby changing the paths of 
trajectories in the vicinity of quasioquilibria. Similarly, control of the not speed of 
motion produced by a model locomotory CPG coupled to an explicit musculoskeletal 
system resulted from the adjustment of CPG trajectories in proximity to unstable 
fixed points of the model [45]. Patterns of biting, swallowing and rejection in Aplysia 
may be understood in terms of sequential traversals between neuromechanical equilib- 
rium points [50, 51]. During normal cell growth and proliferation a living cell passes 
repeatedly through several phases (including cell division), yet the cell "cycle" is typ- 
ically described not as a standard limit cycle but as a sequence of traversals between 
quasiequilibria that act as checkpoints [54].^ In the vicinity of each quasiequilibrium, 
the cell cycle dynamics slows for an indefinite period of time, until a regulating con- 
dition is met [33]. In each of these examples, although the flow strictly speaking 
forms a deterministic limit cycle, the behavior may be actively managed through the 
introduction of variable dwell times that function as control points along trajectories. 

Families of limit cycles verging on a heteroclinic (or homoclinic) cycle such as the 
limit cycles of the one parameter family of systems given by Equations 1.2 provide an 
opportunity for studying the role of saddle points in the control of timing of rhythmic 
behaviors. As a first step, it is natural to consider the structure of the infinitesimal 
phase resetting curves which reflect the sensitivity of the return time along the limit 
cycle to small instantaneous perturbations. When ^ — the flow of Equations 1.2 
does not admit a periodic solution with finite period; consequently the asymptotic 
phase and hence the phase response is not well defined. However, for any > we 
can defined the phase response, and we can study its behavior as the family of limit 
cycles approaches the heteroclinic cycle. 

To fix terminology, for 1 » > 0, consider a trajectory ^{t) following a stable 
limit cycle -f : t £ [0,T) ^ -f{t) € M". Suppose the trajectory is perturbed by a 
small instantaneous displacement, taking x{t~) = j{t + OqT/4) to x{t'^) = ^{t + 
6qT/A) + r. Provided the new initial condition x{t'^) remains within the basin of 
attraction of the limit cycle, the orbit will approach the limit cycle with \\x{t) — 
j{{t + 6\T/A) mod T)]] — >• as t — >■ oo, for some new phase 9\ € [0,4), resulting in 
a shift in asymptotic phase equal to {0\ — 9q) mod 4. The limit cycle's sensitivity 
to weak instantaneous perturbations typically varies both with the phase at which 
the perturbation occurs, 9o{t) = A{{t/T) mod 1), and with the direction in which 
it occurs, T] = f/\\r\\ G M". Differences in sensitivity at different phases, which 
are important for understanding possible control mechanisms, are captured by the 
infinitesimal Phase Response Curve or iPRC, defined as 



As above, is the asymptotic phase associated with a point ^ in the basin of 



^Strictly speciking, the system may be viewed as a limit cycle if the dynamics are embedded 
in a larger space encompassing the control variables as well; the point remains that the periodic 
behavior is strongly influenced by passage near unstable equilibria or near-equilibria, which function 
to regulate the timing of the system. 



Z{0, H, n) = lim 



{9 - <j){j{9T/4:) + eri)) mod 4 



(1.3) 
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attraction for the limit cycle, 9 is the phase at which the instantaneous perturbation is 
applied, /i is the parameter controlling proximity to the heteroclinic, r] is the direction 
of the fast perturbation, 7 is the limit cycle trajectory and T is its period; the latter 
two entities are functions of ^. In order to study the behavior of the phase response as 
the system approaches the heteroclinic configuration, it is important to emphasize the 
limit defining the iPRC in Equation 1.3 is taken before the subsequent limit — ^ 0. 

Analytic calculation of the iPRC is typically accomplished via an adjoint equation 
method, and exact solutions are known in very few cases [23] . In order to perform the 
analysis required to gain qualitative insight into the behavior of the phase response 
under the sequential limits e — )^ 0, — ?• 0, we construct a piecewise linear approx- 
imation to the smooth system. The iris system, described below, is topologically 
equivalent on an open set including the family of limit cycles and the saddle points, 
and qualitatively captures the behavior of the smooth system. For the iris system 
we obtain exact results including the form of the limit cycle for positive /i and an 
explicit formula for the infinitesimal phase response curve. Our main result, stated in 
Theorem 1 below, shows that for the iris system the sensitivity to small displacements 
parallel to the direction of the stable manifold has two regions with distinct sensitiv- 
ity behavior. In the limit, as the family of orbits approaches the SHC, there is an 
interval of phases ip € [0, tpc] for which the infinitesimal phase response goes to zero. 
This interval is separated by a critical phase ipc from a second interval (p G {<fc, 1] 
for which the iPRC diverges to -|-oo. The junction ip = mod 1 corresponds to the 
boundary between regions on which we define a piecewise linear dynamics. We show 
numerically that qualitatively similar results hold for the smooth system. 

Limit cycles in piecewise linear (PWL) dynamical system have been studied pre- 
viously in several contexts. For instance, in the context of Glass networks [20, 25, 26] 
PWL systems have been used to represent the dynamics of idealized genetic regula- 
tory systems. In this case, the structure is somewhat different from that considered 
here, in that the fixed point driving the flow in each piecewise linear region is strictly 
attracting and lies outside the region, rather than having one unstable eigendirection 
and lying inside the corresponding flow region. Consequently, families of limit cy- 
cles verging on a heteroclinic cycle do not appear in Glass networks. PWL planar 
dynamical systems in which a fixed point and a limit cycle coexist do occur in mod- 
els approximating the Fitzhugh-Nagumo equations; such systems were used to study 
traveling wave phenomena [32] and period adding bifurcations under periodic forcing 
[18]. Neither homoclinic nor heteroclinic cycles appear in these systems, however. 
Recently, Coombes investigated phase response curves for limit cycles in both the 
PWL McKean-Nagumo model and a new model related to the Morris-Lecar system 
[16]; this paper exploited the existence of exact solutions for the PRC to study syn- 
chronization in gap-junction coupled networks, in both the strong and weak coupling 
limits. The PWL Morris-Lecar system does contain a homoclinic bifurcation; to the 
best of our knowledge, however, the analysis presented here is the first to obtain exact 
results for the scaling of the infinitesimal phase response curve, as a system of limit 
cycles approaches a heteroclinic or homoclinic orbit. 

2. The piecewise linear iris system. As — )■ 0+, the period of the limit 
cycle in system 1.2 diverges and the asymptotic phase may no longer be defined. 

Nevertheless, the response of the system to small, transient perturbations remains of 
interest. To explore behavior analogous to phase resetting in the /i — > 0+ limit, we 
introduce a piecewise linear planar dynamical system analogous to the smooth system 
in Equations 1.1. Figure 2.1a illustrates the construction: we tile the plane with large 
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(a) (b) 

Fig. 2.1: Construction of the iris system from the smooth system, (a) The flow near 
each saddle point of the smooth system (1.2) is approximately linear in a square region 
aligned with the (orthogonal) eigenvectors. By virtue of the fourfold symmetry of the 
system we can extend each square until it has a side of length 26, remaining centered 
on the saddle. Depending on the extent of the rotational (parametric) perturbation of 
the vector field, there will be square of side a forming a gap around the unstable spiral 
point, (b) We extend the linearized flow for each saddle throughout the corresponding 
square of side 2b, creating a piecewise linear vector field defined on the plane with 
the squares of side a removed. The i*'' square is centered on the i*'' saddle, Xi 
{i = 1, 2, 3, 4). In each square the inward arrows indicate Wi^^ and the outward arrows 
indicate W;"^, parallel to the stable and unstable eigenvector directions, respectively. 



squares of size 26 centered on each saddle and smaller squares of size a between them. 
When the rotation parameter is zero, the large squares align and the square of side a 
vanishes. As the vector field "rotates" in the vicinity of each saddle, the squares rotate 
and a increases from zero. Within each square, we introduce coordinates ^ = (s,u) 
with respect to which the local flow obeys ds/dt = —As and du/dt — u. Here A > is 
the eigenvalue corresponding to the eigenvector (1,0) tangent to the stable manifold 
Wj;^^ = {(s,0)| - 6 < s < 6}. The unstable manifold W = {(0,u)| - 6 < w < 6} 
is tangent to the eigenvector (0,1) corresponding to the unstable eigenvalue, which 
without loss of generality is set to one. The stable and unstable eigenvectors are 
arranged so that when a > the flow on the inner quadrants of each square moves 
clockwise; the flow leaving the square around the i*'' saddle enters the square around 
the i + I'** saddle (mod 4), as shown in Figure 2.1. 

Formally, for k € {1,2,3,4} let the center of the fc*'' square, {xk,yk), be given 
by the real and imaginary parts, respectively, of Zk = \/2(— i)'^(6e^"^/'' + (a/2)e"/^), 
where < a < 6. Define the fc*'' square to be Sk = {{x,y) : \x ~ Xk\ < b &c \y — yk\< 
6}. In the interior of the first square, the flow satisfies x — ^X{x — xi) and y — y — yi. 
The flow in the interior of squares two through four is defined so that the vector field is 
equivariant with respect to fourfold rotation about the origin. Adjacent squares form 
a boundary of length (26— a) . At the boundary between Sk and Sk+i (mod k) , the flow 
leaves Sk and enters Sk+i, by construction. The vector fleld is discontinuous across 
these boundaries. For definiteness, we take the flow at points on the mutual boundary 
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(c) (d) 

Fig. 2.2: The behavior of the iris system depends on the offset a. The system is 
shown rotated through an angle of — tan~^(a/6) for easier visual comparison with the 
smooth system, (a) [a = 0] When the offset a is 0, the unstable manifold of each 
saddle connects to the stable manifold of the next, forming a stable heteroclinic orbit. 
(Thin blue line: sample trajectory.) (b) [a = 0.05] When < a <C 1, a stable limit 
cycle exists that passes close to the saddles. Its basin of attraction extends inward 
to a small, unstable limit cycle around the center (dotted black line), (c) [a — 0.2] 
As a increases, the limit cycle becomes more rounded and the flow along the limit 
cycle more regular, (d) [a = 0.255] If a continues to grow, the stable limit cycle 
is destroyed in a fold bifurcation with the inner unstable limit cycle. Please see 
corresponding movie (file: iris.mpg); this animation cycles through phase portraits 
of the iris system for values of a ranging from to 0.255 and back again. 



Sk n Sk+i of adjacent squares to be defined so that the vector field is continuous when 
approached from square Sk+i, the square into which the trajectories enter. 

The four offset squares may be repeated to form a partial tiling of the plane, 
leaving a complementary set composed of smaller squares of size a. We will view the 
entire system as confined to the 2-torus, however. We leave the fiow unspecified in the 
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interiors of the small squares, and we take any egress points from a largo square into 
a small square to be absorbing, i.e. the flow is set to zero at the boundary of the small 
squares. Since we are interested in the phase response of limit cycles whose basins 
of attraction are entirely contained within the larger squares, the flow in the smaller 
squares is of no consequence. Figure 2.2 illustrates the system, and shows sample 
trajectories for diff'erent values of a. We will refer to the piecewise linear system as 
the iris system because it resembles the iris of a camera, opening and closing as a 
increases or decreases. 

It is clear that the iris system will form a stable heteroclinic orbit when a = 
(corresponding to ^ = in the smooth system); we will show in § 3.2 that it 
also exhibits stable limit cycles passing near the saddles for small values of a > 
(corresponding to /U > 0) as shown in Figure 2.2. When the limit cycle exists, as 
above, we assign a phase 9 to each point 7 in the cycle as the fraction of the cycle's 
period (scaled by 9max = 4) required to reach that point from a defined starting 
point 7(0) = 7o on the cycle, i.e. 

e{i^) = ^mm{{t>0\'y{t) = u}) (2.1) 

where T is the period, i.e. T = min({t > | -f{t) ~ 7o})- It is clear that with this 
transformation d9/dt is constant and equal to 4/r. We define the point where the 
limit cycle first enters the bottom left square (surrounding saddle 6*1) as 70. 

Because the flow is piecewise linear, we may derive analytically the form of the 
phase response curve for any values of A > 1 > a > for which a stable limit cycle 
exists. 

Theorem 1. Let X > 1 > a > 0, let the iris system be defined as in ^2, and 
define the function 

p{u) =u^ -u + a. (2.2) 

1. If the function p has two isolated positive real roots then the iris system has 
a stable limit cycle. Let u denote the sm,allest positive real root of 2. 2. The 
limit cycle trajectory enters each square at local coordinates (1,m) and exits 
each square at local coordinates {s, 1) where s = u^. The period of the limit 
cycle is T = 41og(l/?i). 

2. Let € [0,4) be the phase of an instantaneous perturbation in direction 77 
and let 6 = k + if where (fi G [0, 1) and k G {0, 1,2,3}. If k = then the 
infinitesimal phase response of the limit cycle is 

Z{r], LP, a) = - — -— — — (2.3) 

log(l/u)(u — As) 

where (3{ip) = (^s^^~'^\ u^^ and rj = (77^, 77^) is a unit vector in the Li norm. If 
k € {1, 2, 3} then the infinitesimal phase response is given by Z(r{ , ip, a) where 

r( = T^ri and 7?, = ^ ^ ) ' "^^^ magnitudes of the phase responses to 

perturbations parallel to the stable and unstable eigenvector directions is great- 
est at phases corresponding to boundaries between piecewise linear regions. 

3. As a — > for fixed X, the entry coordinate scales as u = a + o{a), the infinites- 
imal phase response to perturbations parallel to the unstable eigendirection in 
a given square diverges, and the response to perturbations parallel to the stable 
eigendirection diverges when (f & {1 — 1/A, 1); otherwise, it converges to zero. 
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Figure 2.3 illustrates the time course of trajectories of the iris system for A = 2 
and a e {10-3,0.1,0.2,0.24}. 
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Fig. 2.3: Time plots of limit cycle trajectories of the iris system with various values of 
o. When a 4C 1 the trajectories slow dramatically when passing near the saddles and 
travel more quickly between them, resulting in time plots with prolonged dwell times. 
As a increases, the limit cycle moves away from the saddle points and the limit cycle 
trajectory becomes faster and more uniform in speed, resulting in time plots with a 
more sinusoidal character. Note the change in horizontal scale. The horizontal and 
vertical ordinates are shown in black and gray, respectively, with A = 2. From top to 
bottom a = 10-^,0.1,0.2,0.24. Compare Figure 1.2. 



While a general consideration of phase resetting in the vicinity of a heteroclinic 
bifurcation on an invariant circle (or, similarly, near a homoclinic bifurcation) is be- 
yond the scope of this paper, it is natural to conjecture that the limiting behavior of 
phase response curves near such bifurcations may be similar to that observed here. 

Lemmas 3 and 4 in §3 provide necessary and sufficient conditions to guarantee 
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the existence of a stable limit cycle. Lemmas 6, 7 and 8 in §4 develop the response of 
a limit cycle trajectory to a small transient perturbation, leading to direct calculation 
of the infinitesimal phase response curve. Lemmas 10 and 11 in §5 describe the 
asymptotic behavior of the iris system in the heteroclinic limit, i.e. as a — )■ 0. In §5 
wc also compare analytic and numerical results for the phase response curves of the 
iris system. In §6 we numerically explore the isochrons of the iris system. Finally 
in §7 we numerically obtain phase response curves for the smooth system given by 
Equations 1.1 and compare their structure with those of the iris system. 

In [11] Brown et al. studied phase response curves for limit cycles near the four 
codimension one bifurcations leading to periodic firing in standard neuronal models 
(saddle-node bifurcation of fixed points on a periodic orbit; supercritical Hopf bifur- 
cation; saddle- node bifurcation of limit cycles; homoclinic bifurcation). Their analysis 
of the homoclinic bifurcation corresponds to our analysis of the iris system in a certain 
limit; for a detailed comparison see §8.2. 

3. Limit Cycles in the Iris System. We now prove several results about the 
iris system that we will need to prove Theorem 1. We start by studying the trajectory 
within one of the square regions to construct a map from the time and position of entry 

into the region to the time and position of egress out of the region. Next, we connect 
four of the linearized regions together. We then prove the existence of a limit cycle 
for sufficiently small values of o > 0. We also examine the effects of a perturbation of 
the trajectory within the neighborhood on the exit time and exit position. We use the 
maps thus derived to show that a stable limit cycle exists, and use the perturbation 
results to determine the asymptotic effect of a small perturbation on the phase of the 
oscillator [i.e. the phase response curve). First, we will examine the dynamics within 
a single linear region around a saddle of the iris system. 

Remark 2 (Nondimensionalization) . Assume each saddle of the iris system has 
two real orthogonal eigenvectors, one stable and one unstable. Assum,e the region 
around each saddle is a square aligned with these vectors and centered on the saddle 
with a length of2L along each side (so that the saddle is a distance of L away from 
each edge). We will refer to the unstable eigenvalue of the saddle as A„ > and, the 
stable eigenvalue as Xg < 0. The trajectories within the region have the dynamics 

ds du /r, i\ 

- = - = Xuu. (3.1) 

Assume that when a trajectory leaves the edge of one region atpositionxf = (s/,u/) = 
{sf,L), it enters the next region with an offset of a, i.e. Xi = {si,Ui) = {L,Sf + a). 
We assume that there is a cycle of four regions the trajectory may traverse this way. 

We can nondimensionalize the system by defining the new state variables s* = s/L 
and u* = u/L, by rescaling time as t* = Xut and the offset as a* = a/L. We may 
define X = — As/A„ as the saddle ratio, which will play a critical role in stability (^3.2). 
Using these definitions, the governing equations become 

— = -As*, = u*. (3.2) 

dt* dt* ^ ' 

We will use the non-dimensionalized version of the system for the remainder of the 
paper. For simplicity of notation, we will omit the asterisks in the sequel. 

3.1. Dynamics Within A Linear Region. Within a single square region sur- 
rounding a given saddle point, we will use the local coordinates {s,u) to represent 
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the displacement parallel to the stable and unstable eigenvectors, respectively. The 
line segments {(s,0) G M^lO < s < 1} and {(0,u) € M^lO < u < 1} arc part of 
the stable and unstable manifolds of the saddle, respectively, forming separatrices of 
the flow in the square. We will restrict attention to flow entering the square along 
the edge € M^|0 < u < 1}. It is easy to see that a trajectory entering at 

an initial point Xi = (si,Uj) = at time ti with Ui > will exit at point 

{sj.uj) = = at time log(l/Mi). We deflne the map /j : M — >■ M from 

an entry position along the ,s = 1 edge to an exit position along the u = 1 edge of 
the region of linear flow, and a function Ti{ui) representing the transit time from an 
entry position 

fi{ui) = Ti{ui) = log(l/u,)- (3.3) 

Note that /; is a monotonically increasing function of Wj and that Ti is a monotonically 
decreasing function of Wj. 

The closest approach of a trajectory to the saddle point occurs when the position 
X = {s,u) is perpendicular to the velocity v = {—Xs,u), or = — Aexp[— 2At] + 
uf exp[2t\. Similarly, one may calculate the time at which the speed of the trajectory 
is minimal. These times are 

log A — 2 log Uj 3 log A — 2 log Uj 

^closest = 2(A + 1) ' ^slowest = 2(A + 1) ' ^ 

We can examine the position of closest approach by integrating the system from the 
entry position for a time i^iosest ' which gives the location 

"closest = Wje*<^'°='==' = u, (A/u^)^^^^^"'"^^ 
Taking the ratio of the two coordinates we find 

"closest = A ^ ^closest) (3*5) 

and thus the point of closest approach of each trajectory in this quadrant lies along 
a line passing through the saddle. A similar calculation using ^slowest shows that the 
slowest points lie along the line 

"slowest X ^ ^^slowcst (3-6) 

which also passes through the saddle point at the origin of the local coordinate system. 

3.2. Dynamics Across Regions. We will next explore the conditions under 
which the iris system contains a limit cycle. 

Lemma 3. The iris system described in ^2 contains a limit cycle iff the function 

p{u) = — u + a 

has isolated real roots G (0, 1), with <u^. These roots exist iff X> 1 and 

AV(i-A)_AV(i-A)+a<0. 

Proof. To track crossings between regions, we will add a subscript to the variable 
names indicating how many region crossings have occurred since time t = 0, so that 
Ui^s is the ingress location (uj) immediately after the third crossing. 
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When a > 0, trajectories leaving one region enter the next with an offset o, 

i.e. (si^n+ijUi^n+i) = (Ij -S/.n + Q^)- Wc therefore have a second monotonically increas- 
ing map from the point of egress afong the edge of one square to the entry point along 
the edge of the next square, /e : sy_„ — >• Ui^n+i = + a- 

We can now form a map from the entry position along the edge of one linear 
region to the entry position along the edge of the next region via the composition 

h{ui) = (/e o fi){ui) = (Ui)^ + a. (3.7) 

As the composition of two monotonically increasing maps, h is also monotonically 
increasing. Noting that the entry edges form a transversal section of the flow, it is clear 
that h^, the fourfold composition of h, forms a Poincare map for any cycles crossing 
this edge. Limit cycles will form isolated fixed points in this map, so to find potential 
limit cycles we look for points where Ui = h'^iui). Because h is composed of four 
identical monotonic maps, the limit cycles will also be fixed points in these component 
maps and thus we have a fixed point at Ui = u if and only if u = h{u) = + a, or 
equivalently p{u) = — u + a = {cf. Equation 2.2). 

We now consider the potential values of A, which by definition must be positive. 
First, consider < A < 1. Rewriting p(u) as u^{l — u^~'^) + a and remembering 
that u S (0,1), we can see that p{u) is now the sum of a positive product and a 
non-negative parameter, and thus cannot be zero. Next, consider A = 1. In this case 
p{u) becomes just a, and thus p = implies a = 0. Note that p is independent of u 
and thus every u G (0, 1) is a solution, corresponding to a one parameter family of 
neutrally stable, non-isolated orbits. Therefore A must be greater than one if a limit 
cycle exists. 

Differentiating p twice with respect to u, we find that p is convex for all u > 
when A > 1, and differentiating it once we find that p has a minimum at 

(uiUn = AV(i-A). (3.8) 

Because p{u) is continuous, and equal to a when u = or m = 1, by the midpoint 
theorem it will have a root between and (ui)min and a root between {ui)inin and 1 
iff 

> r{{uiUn) = («i)min ' {^iUn + a = X^/(^-^^ - aV(i-^) + a. (3.9) 

We will call the smaller root and the larger (possibly identical) root u^. □ 
Figure 3.1 illustrates the map h = f^o fi when A = 2 and a = 0.2. 
We now examine the stability of the roots and u*. 

Lemma 4. If the roots v) and of Equation 2.2 exist and are distinct, gives 
the entry position of a stable limit cycle along the s = 1 edge of a region and gives 
the entry position of an unstable limit cycle along the same edge. 

Proof. We can determine the stability of the two points by examining the deriva- 
tive of h at these points; if \dh/dui\ > 1 the fixed point will be unstable, and if 
\dh/dui\ < 1 the point will be stable. Because the existence of the roots implies 
A > 1, we know dh/dui ~ {ui)^~^\ will be positive. Differentiating h a second time 
we see that dPh/dul = {ui)'^~^{\ — 1)A is always greater than zero, and thus dh/dui 
is monotonically increasing. Since a fixed point in a map is stable iff the magnitude 
of the derivative at that point is greater than one, we find the point where dh/dui 
passes through 1: 

1 = dh/dui = A(ui)^"^ 
13 
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Fig. 3.1: The map h — fe o fi from the entry position of a trajectory along the edge 
a square to its entry position along the edge the next square forms a map analogous 
to a return map (blue curve). This map is continuous and monotonically increasing, 
so a limit cycle can only exist where a trajectory enters the next square at the same 
position it entered the current square (green line indicates the identity map). In this 
example there are two intersections corresponding to the stable and unstable limit 
cycles. Here a — 0.2 and A = 2. The blue curve meets the vertical axis at the offset 
between squares a, and changing the offset raises or lowers the blue curve without 
changing its shape. Varying a allows for 0, 2, or 1 limit cycles (corresponding to 
Figure 2. 2d, Figure 2.2c, and the fold bifurcation that occurs between them). 



or 

This, however, is just (Mi)rnin from Equation 3.8. Because this lies between the two 
roots, dh/duil^^f < dh/ dui\^^^:^ . = 1, and thus is a stable fixed point of /i, and 

1 = d/i/c?Ui|(jj.-)^ < dh/dui\^t and thus is an unstable fixed point of h. Because 
/i is a Poincare map at the entry plane of the region, these fixed points correspond to 
the entry points of a stable and unstable limit cycle respectively. □ 

Remark 5. Lemmas 3 and 4 together establish Theorem 1, Part 1. 

As just shown, for small positive a two limit cycles exist - one stable and one 
unstable. As a — >■ 0, the stable limit cycle is destroyed in a heteroclinic bifucation 
and the unstable limit cycle collapses to a point at the center of the system. As a 
increases, the two limit cycles collide in a cycle-cycle fold bifurcation. The collision 
occurs when the inequality (3.9) becomes an equality. The two roots converge to 
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trajectories spiral into center 




A 



Fig. 3.2: Two-parameter bifurcation diagram for the iris system. When the stable- 
to-unstable eigenvalue ratio A is sufficiently large relative to a, a stable and an un- 
stable limit cycle coexist (region labeled 'stable limit cycle'). As o increases the two 
limit cycles are destroyed in a fold bifurcation (line labeled 'fold bifurcation'). As a 
approaches 0, the system approaches a heteroclinic orbit (line labeled 'heteroclinic 
orbit'). When a is sufficiently large or A < 1, trajectories spiral into the center as 
shown in Figure 2. 2d (area labeled 'trajectories spirals into center'). When A = 1 and 
a = 0, the system becomes a square filled with neutrally stable orbits (point labeled 
'neutrally stable orbits'). 



(W'i)min sls thc two limit cycles that cross the section at the two roots merge. Figure 
3.2 shows thc bifurcation diagram. 

4. Effects of a small instantaneous perturbation. We now examine the 
effects of a small perturbation of a trajectory starting on the limit cycle. We will 

first explore thc effects on the exit time and exit position from the square, and then 
extend our analysis to the asymptotic effect of the perturbation on phase over many 
cycles (the phase response curve or PRC). 

4.1. Initial effects of a small perturbation. Assume a trajectory begins on 
thc stable limit cycle at position Xq = (1,m^) at time t = 0, and receives a small 
perturbation of size r in the direction of a unit vector rj — {rjmrja) at time t' . We 
will consider only perturbations that occur before the trajectory leaves the square 
{i.e. < t' < ri(ut) = tI) : earlier or later perturbations can be reduced to this 
case by an appropriate time shift and discrete rotation. We will also require that 
the perturbation not push the trajectory out of the basin of attraction of the limit 
cycle. The infinitesimal PRC we calculate below is obtained in the limit of small 
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perturbations; for finite a > 0, the basin of attraction has finite width, and trajectories 
that escape the limit cycle's attracting set arc not assigned an asymptotic phase. 

It is straightforward to calculate the effects of the perturbation on the time and 
position at which the trajectory leaves the square. We first consider the typical case 
where the perturbation does not push the trajectory across the border of the square. 

Lemma 6. Consider a trajectory initially on the stable limit cycle 7 of the iris 
system such that 7(0) = {l,u^), and a perturbation Ax of size \\Ax\\ = r in the 
direction of unit vector 77 = {r]s,riu) = Ax/r at time < t' < tI . Assume the 
perturbation does not push the trajectory out of the basin of attraction or into another 
square. Then the perturbation will result in a change in the position of entry to the 
next square 

Au[ = (r?,(«te*y + r?„A(wt)A-ie-*')r + o(r) (4.1) 

and a change in transit time of the square 

AT[ = -{r,ue-*'/u^)r + o{r), (4.2) 

as r ^ 0. 

Proof Immediately before the perturbation occurs, the trajectory will be at the 
point 

lim x(t) = (e-^*',wtet'). 

The perturbation shifts the position to 

lim x{t) = {e~^*' + rr]s,uh*' + rrju). 
t->-t'+ 

The trajectory will exit the square once x^ grows to 1 at time 

T[^t'+ (- log(u^e*' + r7?„)) 

and enter the next square along the s = 1 edge at position 

u[=a+ (e"^*' + rr?s)e-^('^i'-*'' =a+ (e"^*' + rT]s){uh*' + rr?„)^ 

For r <C 1, we may write 

Ti=t' - \og{u^e*') - {r]ue-'' /u^)r + o(r), 

«; =a+(u^)^ + (r7,(w^e*y + r;„A(u^)^-ie-*')r + o(r) 

as r — >■ 0. Noting that t' — log(w'''e* ) = — log(ui) = tI and a + {u'^)^ = , these 
expressions may be viewed as the exit time and position for the limit cycle with a 
linear correction and additional higher order terms in r. That is, 

Ti=Tl-{r]ue-''/u^)r + o{r), 

u[ =u^ + (%(Mte*y + r?„A(ut)^-^e-*')r + o(r). 

Subtacting tI and respectively then gives our result. □ 

Perturbations that push the trajectory across the boundary between two squares 
can be handled without much additional difficulty. For example, we can treat a 
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perturbation that crosses into the next square as an advancement of the trajectory 
to the point that it enters the next square followed by a perturbation to the new 
location in the new square. Perturbations to the previous square can be handled in a 
similar fashion, as can perturbations to the diagonal square. However, for any finite 
a > and any perturbation time other than t' — or t' — tI, the perturbed point 
will remain in the same square as the unperturbed point for sufficiently small r > 0. 
Therefore the additional cases do not alter the calculation of the PRC except possibly 
at a set of finite points, and they will be omitted. 

4.2. Subsequent effects of a perturbation. In the absence of perturbation, 
the limit cycle passing through 7(0) = will exhibit a (constant) sequence of 

edge crossing locations = u^,n = {0, 1,2, •••} at a sequence of crossing times 

= uTI at constant intervals AT^ = T|. A single perturbation Ax = r{r]s,ilu) at 
time t' will lead to a new sequence of edge crossings {u'^ „} and crossing times {T^}. 
Following the perturbation, the time spent between the n*'' and n + 1"* crossings is 
ofi^set from the unperturbed dwell time: T^^-^ —T^ = tI + AT^^^. Next we calculate 
the change in entry position for subseqTicnt edge crossings (Au^ = u'^ ^ — u^) and the 
dwell time offsets for subsequent squares {AT^_^_i). 

Lemma 7. Under the conditions of Lemma 6, the entry position after the n*^ 
crossing [n > 1) will he offset by 

A< = (A(ut)^-i)"-' (%(w^e*y + r?„A(nt)^-ie-*> + o(r), (4.3) 

and the time spent in that square will be offset by 

AT^+i = (t.t)-i (A(ut)A-i)"-i (%(ute*V + r,uX{u^)^-'e-'')r + o{r), (4.4) 

as r ^ 0. 

Proof. A perturbed trajectory will enter the r?,"' edge with an offset Au'^. This 
situation is equivalent to a trajectory entering along the limit cycle that experiences 
a pcrtm'bation in the u direction (r/ = (0, 1)) of magnitude r = Au'^ immediately 
upon entering the square. We can thus use Equation 4.1 with t' = to find the entry 
position along the next (n + 1"*) edge: 

Au'^^,=X{u^)'-'Au'^ + o{Au'J. 

Therefore to leading order the sequence of offsets follows a geometric series, with 

closed form 

A< = (A(M^)^-i)""' Au[ + o{Au[) 

= (A(nt)^-i)"-^ (%(«^e*y + r?„A(wt)^-ie-*')r + o(r), 

as r — > 0. Substitution into Equation 4.2 provides the deviation in the time of passage 
through each square, compared with that for the limit cycle trajectory. Again setting 
T] = (0, 1), r = Au'^, and t' = gives 

Ar;+i = -(z.t)-iA< + o«) 

= -(wt)-i (A(i,t)A-i)"-i (r?,(ute*y + r,„A(«t)A-ie-*> + o(r). 

□ 
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We now examine the asymptotic eflFect of the perturbation at long times. Because 
u' is less than one, both Aw, and At will approach zero for large n. This is equivalent 
to saying that the orbit will asymptotically return the stable limit cycle, as expected. 
It may, however, return with a different phase than the unperturbed trajectory. 

Lemma 8. Under the conditions of Lemma 6, the cumulative change in crossing 
times after n crossings is 

^ .t-A(.t)A ^ + oir) 

(4.5) 

which converges to 

AC =-^^M^!^l^±^^r + oM (4 6) 

as n — > oo. 

Proof. The cumulative change in crossing times is the sum of the dwell time 
offsets in the previous squares: 

n n 

AC'„ = J2^U = ATi + Y,AU 

fc=l k=2 

Substituting in Equations 4.2 and 4.4 yields the result: 



-1 



AC; = -'^r + {X{u^)'-')'~" iVsiu^e*')' + ^„A(«t)A-ie-*') + o{r) 

fe=2 

Vue-'' + (r?.(«te*y + r?„A(t.t)^-ie-*') (A(ut)A-i)^'j ^ ^ ^^^^ 
T [vue-^' + iVsiuU^Y + r/„A(t.t)A-ie-*') ^^^^^y^^) r + o(r) 

r + o(r). 



i/,t - A(ut)A 

Recalling that < < 1 and A > 1, taking the limit as rn- oo gives 

«t-A(«t)A ^ + 

as required. □ 

The limit cycle enters each square at location (1, u^) and exits at ((u^)^, 1), so it 
is natural to define the exit coordinate along the stable eigenvector axis as s^ = (u^)^. 
Recalling that = + a, Equation 4.6 may be simplified as 

, r?s(sV/*') +77„e-*' r]-l3 
AC^ = r + o{r) = ~^^^r + air). (4.7) 

The vector /3 = ^s'^e^* , e~* ^ may be thought of as a solution to the same differential 
equation as a trajectory within the first square, ^ = (s,u), but with time reversed 
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and initial conditions set to the far end of the square. That is, 



/-AO 



fx \ ^^^^ f .t 



dt V -1 / ' V 1 

4.3. Infinitesimal phase response curve. Wc now derive the infinitesimal 
phase response curve. Initially, we restrict attention to perturbations occurring within 
the first quarter cycle, i. e. within the first square. The PRC for perturbations occur- 
ring after k additional border crossings is obtained via a discrete rotation operation. 

We define the phase of a point on the limit cycle as in Equation 2.f, i.e. as 
the smallest amount of time required to reach that point from the point where the 
limit cycle enters the first square. Therefore the asymptotic shift in crossing times 
translates into an asymptotic shift in phase as 

4 ^ ^, AC^ 
Aq3 = —AC = — — . 

^ T log(l/^it)' 

recalling that T = 4T^ = 41og(l/u'''). Since e~* = (u^)'^ we may write /? in the form 
= = ((nt)^(i-*'), {u^r) = ((st)a-*'), («t)c.) , (4.8) 

simplifying the subsequent analysis. Thus the asymptotic change in phase for a per- 
turbation at phase (p G [0, 1) is 

^^= log(l/IiKlt-A.t) ^ + ^(^)- 

The infinitesimal phase response curve, obtained by taking the limit of A(p/r as r — )• 0, 
is therefore 



log(l/wt)(ut-Ast)' 

which is Equation 2.3. 

For a perturbation at phase lys > 1, let be the number of boundary crossings 
preceding the perturbation, i.e. the positive integer satisfying ip & [k,k + 1). Define 
a fc-fold quarter rotation operation on the phase variable and the direction of the 
perturbation by 

(fi ^ (fi' ^ ip- k (4.9) 
r]^r]' = n'^r] (4.10) 

where TZ = ( ^ n ]■ Then the infinitesimal phase response due to a perturbation 



1 ^ 

in direction 77 at phase 6* > is given by Z{ri' , ip' ,a). 

The dependence of the PRC on the phase is entirely contained in the term 7] ■ 0. 
Consider the components of the response to perturbation in the horizontal and vertical 
directions. When tp € [0, 1) the component 0s along the direction of the stable 
eigenvector corresponds to Px, the component in the horizontal direction. Similarly 
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the components along the unstable direction ^„ and the vertical component Py are 
identical. Moving clockwise around one full orbit we observe that: 



when </?e[0,l), 13^ = 13, 

when ip G [1,2), ji^ = /3i 

when ip e [2,3), (3x = 

and when ip G [3,4), = — 



'u 



s and py = Pu\ 

u and Py = -/3s; 

/3s and /3y = 

and Py = Pa- 



For each component the resulting phase response curve consists of two continuous 
segments separated by a jump of size {u^ + s^), with antisymmetry /3((/? + 2) = — /3((p). 
The horizontal response component, P^, peaks at Pxi'^) = 1 and is strictly positive for 
e [0, 2). Similarly, (iy peaks at /3j^(0) = 1 and is strictly positive for ip e [0, 1)U [3, 4). 
Hence the extrema of both components occur at phase values corresponding to the 
boundaries between adjacent squares. 

Remark 9. This concludes the proof of Theorem 1, Part 2. 

Figure 4.1 shows examples of the full iPRCs as a function of ^ € [0, 4) for A = 2 
and a e {10-3,0.1,0.2,0.24}. 

5. Asymptotic phase resetting behavior as a ^ 0. We wish to study the 
asymptotic behavior of the phase resetting curve at the heteroclinic limit, that is, 
as the parameter a ^ 0. This limit corresponds, at least by analogy, with the limit 
/i — > in the smooth vector field system of Equations 1.2. 

Lemma 10. Under the conditions of Lemma 6, the limit cycle entry value scales 
as = a + o{a), as a^Q. 

Proof. Recall that ti) is the smaller root of /o(u), so u''' = a + {u^)^. Since A > 1, 
clearly = a + o{v)), as v) — >■ 0, and > a whenever > 0. For sufficiently 
small a > we have as an implicitly defined function of a. It is straightforward 
to sec that du' /da = 1 when = 0. We wish to show that u'' ^ a + o(a) as well, 
i.e. that lima^o{u^ — a)/a = 0. Since = a + o(a), for any e > we can find 
a ^ > such that 1 — a/v) < min[e/2, 1/2] whenever < < ^. If a < ^ then 
{v) -a)/a = 1/(1- (l-a/ul')) - 1 = 1/(1 - .x) - 1 = x/{l~x) where x = 1-a/ut. 
Since x < min[e/2, 1/2], we are guaranteed that x/{l — x) < e, as required. □ 

As a decreases, the peaks of the iPRC grow, while their width shrinks. The 
integral under each strictly positive region is 



which goes to zero as — )• or equivalently as a ^ 0. Lemma 11 characterizes the 
behavior of the "normalized" iPRC a{(f) = f3{(p)/V as a 0. 

Lemma 11. For the iris system the vector a{(p) = P{ip)/V has the following 
properties as a ^ 0: 

1. For (f = 1 or 3, ax ^ ±oo, respectively. For (p ^ {1, 3}, ax — > 0- 

2. For ip = or 2, ay ^ ±00, respectively. For tp ^ {0, 2}, ay — >■ 0. 

3. /q ||a((/?)|| dip = A (in the Li norm). 

Proof As a -> both V ^ and u'^ 0. At (p =^ 1, Px = 1 for all a > 0, 
so as a — >■ the ratio I3x{l)/V +00. On the other hand, for fixed ip e [0, 2)\{1}, 
(^yt)A(i-(,3) _s. as — > 0. The behavior of ax or ay within the other intervals follows 
from the symmetry relations discussed in the preceding paragraph. Thus 1 and 2 are 
established. For 3, note that by construction each interval of length two centered on 




(5.1) 
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Fig. 4.1: Infinitesimal phase response curves (iPRCs) for the iris system with vari- 
ous values of a. As a becomes small and the limit cycle approaches the heteroclinic 
orbit, the phase response curve becomes dominated by peaks at the edges between 
squares where flow changes from compressing trajectories outward to expanding them 
inward. As a grows and the flow along the limit cycle becomes more uniform, the 
iPRC becomes less sharply peaked. Line (dark blue): analytically obtained iPRC 
for a perturbation in the positive x direction. Points (dark blue): numerical iPRC, 
calculated using an instantaneous perturbation of 10"'' in the horizontal direction. 
From top to bottom, a = 10"*^, 0.1, 0.2, 0.24. The light gray curve represents the 
corresponding iPRC for perturbations in the vertical (positive y) direction. The dis- 
continuities between the positive and negative portion of each curve are steps of height 
{u^ + s^) (see §4.3 for further details). 



a peak of one component makes a unit contribution to the integral. That is, 



1= / axi(p)dip = - / axi(p)d(p^- / ay{ip) dip = / ay{ip)d(p. 

J[0,2) J[2,4) J[l,3) i[3,4)u[0,l) 

The integral of the Li norm ||q;(0)|| is the sum of these terms. □ 

Remark 12. Lemma 11 and the symmetry of a together imply that each compo- 
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nent of a converges weakly to a sum of delta function distributions on the circle: 



a^i'f) ^ d{f - 1) - - 3) (5.2) 
a,(v5)^%)-%-2), (5.3) 

as a — >■ 0. 

Rewriting the full PRC in terms of a gives 

^' = log(V.?(:i!A.t) = • -(^)) ^("^) (5-4) 

where M{u'^) = )^(^i^g^/'~f^yi ' '^^^ function M represents the magnitude of 

the PRC and diverges as (m'^ log^(u''')) ^ as — >■ 0. 

The asymptotic phase response to an arbitrary instantaneous perturbation in 
direction 77 is a sum of the response to the perturbation component along the stable 
eigenvector direction, t]s and along the unstable eigenvector direction, rju, within 
whichever square the perturbation occurs. The asymptotic behavior of the phase 
response in the heteroclinic limit for each component is distinct. We rewrite the PRC 
to emphasize these contributions thus (restricted, for convenience, to the SW square, 
or [0,1)): 

^^"''"'"^^ log(l/nt)(.t-Aat) • ^'-'^ 

First consider the response to perturbation along the unstable direction. Recall 
that for any p > 0, u^logu — >■ as u — ^ 0+. It follows after a brief calculation that 
for any (p € [0, 1) 



lim 



log(l/'u)('U — Xu^) 



= +00. 



This result means that regardless of the phase of the perturbation, the sensitivity 
of the asymptotic phase to small displacements parallel to the imstable eigenvector 
diverges as the system approaches the heteroclinic bifurcation. The result is intuitively 
appealing because when the system is close to the heteroclinic bifurcation the long 
period leads to a larger "compounding" effect enhancing the cumulative result of a 
small perturbation within a given square. However, as is often the case when dealing 
with nested limits, intuition can be misleading. In contrast to the preceding situation, 
consider the response to perturbation along the stable direction. The behavior of the 
limit now depends on the parameter A: 



lim 

U-¥0+ 



log(l/u)(u — Xu^) 



0, e [0, 1 - i/A] 
+00, (/? e (1 - 1/A, 1) 



Thus, in the heteroclinic limit, small perturbations parallel to the stable eigenvector 
become inconsequential if they occur at an early enough phase. At a particular phase, 

<^* = 1-^, (5.6) 

the response becomes highly sensitive to arbitrarily small perturbation. 
Remark 13. This concludes the proof of Theorem 1, Part 3. 
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For comparison we can obtain from Equations 3.4 the fractional phases at which 
the trajectory is closest to the saddle point (in the first square) or moving at the 
slowest speed: 

_ l-logA/(21ogut) _ l-31ogA/(21ogMt) 

'^closest ~ X + 1 ' '''slowest ~ A~+l ' 

As u''' — >■ both phases converge to fo = 1/(A + 1), which is distinct from the phase 
at which the PRC's asymptotic sensitivity to perturbations in the stable direction 
changes. Curiously, the asymptotic value ifo of the phases of the slowest and closest 
point along the trajectory coincides with the critical phase for sensitivity, precisely 
when A is equal to the golden ratio, (1 + •\/5)/2 « 1.618- • • , because under those 
conditions 1 - 1/A = 1/(A + 1). 

6. Isochrons. Provided a > in the iris system, we may define isochrons as 

the level sets of the asymptotic phase function 0{x), for any point x in the basin of 
attraction of the limit cycle. As described in [11, 28, 31] and elsewhere, the points 
on a given isochron will converge over time to a single point on the limit cycle with 
a particular phase. Thus, in a sense, they represent points with the same asymptotic 
"time". As discussed above, given a limit cycle {j{t)}JlQ the asymptotic phase of a 
point xo is defined as the unique value in ^ e [0, ^max) such that the limit 

lim |x(0-7(i + %o))| = 0, (6.1) 

where x{0) = xq is the initial condition for the trajectory. The usual approach for 
finding isochrons is based on an adjoint equation method, cf. ([11], Appendix A; see 
also [23], Chapter 11). Prom the chain rule, the phase field must satisfy 



.,XW,= [(W),XW)].| = ^ (6.2) 



T{a) 



where T(a) is the period of the limit cycle for a given value of a > 0. Let (0 < s < 
1,0 < M < 1) be local coordinates relative to any saddle of the iris system, i.e. such 
that the flow satisfies s = —As, u = u. Writing 9^ for 80/ dv, Equation 6.2 reads 

-As^, + uOu = ema^/T{a). (6.3) 

By virtue of the fourfold rotational symmetry of the iris system, the following bound- 
ary conditions are required in addition to the PDE 6.2: for < s < 1 — a, 

e{s,l) = 6{l,s + a) + l (6.4) 



The boundary condition renders the PDE nonlinear, and a simple closed form solution 
is not available. Solving the PDE numerically, however, is equivalent to finding the 
isochrons via direct simulation of trajectories for a suitable grid of initial conditions. 
We used a 400x400 square mesh of initial conditions within each square, tracking 
solutions until either they left the system of large squares (initial conditions on the 
interior of the unstable limit cycle) or converged with a suitable tolerance to a small 
neighborhood of the limit cycle, at which point the asymptotic phase could be as- 
signed. Figure 6.1 shows an example of the isochrons generated with this technique. 
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Fig. 6.1: Isochrons for the iris system obtained via numerical iteration of the map 
from initial conditions to a sequence of boundary crossing. The iris system on the 
torus has two distinct stable limit cycles for a > 0, illustrated here in two color 
schemes (green-blue and fuschia- violet). The tiling of the plane by periodic extension 
of the basic torus system is shown. Solid lines indicate the stable limit cycles, dashed 
lines indicate the unstable limit cycles. All limit cycles rotate clockwise. The flow 
inside the small squares is not defined, and points inside the unstable limit cycles are 
not part of the basin of attraction of the stable limit cycles. The stable manifolds 
of each saddle comprise the separatrices between the basins of attraction of distinct 
stable limit cycles. The system shown has a — .2 and A = 2. Colors (online) indicate 
the phase. Note the compression of the level curves near the boundaries between 
adjacent squares, which correspond to peaks in the phase response curve. Note also 
that the isochrons are not strictly parallel to the stable eigenvalue direction (c./. the 
discussion in §8.2). Please see corresponding movie (file: iris_isochrons .mpg); this 
animation illustrates the flow together with the isochron structure, by evolving points 
of equal phase forward together in time. After five seconds the animation changes to 
showing isochron structures for different values of the unstable/stable manifold offset 
parameter a, from a = to a = 0.247 and back again. 
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7. Smooth System. In order to compare the family of infinitesimal phase re- 
sponse curves obtained for the iris system with that of a family of continuous vector 
fields experiencing a similar heteroclinic bifurcation with fourfold symmetry, we nu- 
merically evaluated the phase response for small instantaneous perturbations in the 
horizontal and vertical directions for limit cycles given by the system 1.1 for positive 
values of the twist parameter ji. Figure 7.1 shows good qualitative agreement with 
the concentration of phase response sensitivity at points along the trajectory near 
integer values of the phase ip, cf. Figure 4.1. 

8. Discussion. The iris system introduced here is one of only a handful of non- 
linear dynamical systems for which an explicit form for the limit cycle or the phase 
response curve is known. For this system we show analytically that the sensitivity of 
the response to small instantaneous perturbations depends (as one would expect) on 
the distance to the heteroclinic bifurcation, on the phase at which the perturbation 
occurs, and the direction of the perturbation in the phase plane. In addition, we find 
the intuitively appealing result that regardless of phase, as the static perturbation 
away from the heteroclinic vector field diminishes, the phase response as measured by 
the iPRC becomes hypersensitive to any small perturbation parallel to the unstable 
eigenvector direction, in that the iPRC diverges as the bifiircation parameter of the 
iris system a — )■ 0. Unexpectedly, however, we also find that for perturbations along 
the stable eigenvector direction, the response to some will diverge while the response 
to others will have no effect in the limit as a — > 0. The system appears ultrasensitive 
in this case to the phase at which the perturbation occurs, with perturbations suffi- 
ciently far in advance of the approach towards the saddle point effectively absorbed 
by the flow, in contrast to perturbations beyond a critical phase, ip.^ = 1 — 1/A, which 
diverge as a — )■ 0. From a biological point of view, this dual sensitivity to the timing 
and direction of the transient perturbation is significant if flows structured by flxed 
points are to be exploited in nature as control points for rhythmic behaviors. In neiiral 
systems, for instance, it is commonplace to consider perturbations restricted to a sin- 
gle dimension in a multidimensional flow, namely perturbations along the dimension 
of membrane potential [11]. When a neuron's membrane potential lies in the linear 
regime, its dynamics is typically dissipative due to membrane conductances, tending 
to align the voltage direction with the stable direction for a subthreshold fixed point, 
at least for some portion of the flow. In general, if slowly adaptive modulatory pro- 
cesses within an organism's control circuitry can impose quasistatic perturbations of 
an existing vector field to move trajectories closer to or farther away from a homo- 
or heteroclinic bifurcation point, the modulatory processes can filter which perturba- 
tions the system becomes sensitive to at different phases of the ongoing oscillation. 
Investigating these possibilities in specific systems such as the feeding central pattern 
generator of the marine moUusk Aplysia will make it possible to test this prediction 
empirically. 

8.1. Sensitivity and control. Our results suggests that for systems in which a 
control parameter can change the approach of trajectories towards a hyperbolic saddle 
fixed point with a one dimensional unstable manifold, the sensitivity to specific kinds 
of perturbations could be actively managed by manipulating the critical phase ip^. 
This is a novel kind of control for rhythmic behaviors; additional mechanisms for 
control include manipiilating the period of a typical trajectory by controlling the 
time spent near a saddle point, and regulating the variable dwell time in different 
states for a limit cycle trajectory passing near multiple unstable fixed points. 

Oscillations arising from co-dimension one bifurcations other than a saddle node- 
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Fig. 7.1: Infinitesimal phase response curves for the smooth system with various values 
of /i. As /i becomes small and the limit cycle approaches the heteroclinic orbit, the 
phase response curve becomes dominated by peaks much like the iris system shown in 
Figure 4.1. As grows, the phase response curve becomes a sinusoidal as the system 
approaches an Andronov-Hopf bifurcation. Points connected by dashed line (dark 
blue) : numerical infinitesimal phase response curve, calculated using an instantaneous 
perturbation of 10~^ in the horizontal direction. Dots (dark blue): infinitesimal phase 
response curve for perturbations in the horizontal direction. Light gray: infinitesimal 
phase response curve for perturbations in the vertical direction. From top to bottom, 
fi = 10-3,0.1,0.3,0.45. 



homoclinic have phase response curves near the bifurcation point that are smooth 
and sinusoidal. By contrast, many phase response curves observed in real systems 
show sharp transitions between small and large responses, even in systems generating 
highly regular behaviors such as swimming in the lamprey. Lamprey motor units 
participating in the CPG underlying smooth swimming showed broad regions of near 
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zero phase response combined with steep changes in sensitivity in other regions, going 
from zero to peak phase response in an interval corresponding to about 10% of the 
limit cycle duration (c./. Figure 6 of [55]). Preliminary data obtained from the marine 
mollusk Aplysia californica during feeding behavior show extended, variable dwell 
times in preferred regions of phase space separated by rapid transitions from region to 
region [44], suggesting the possibility that this central pattern generator's dynamics 
may resemble that of a limit cycle making close encounters with multiple unstable 
fixed points. 

The infinitesimal PRC is equal to the gradient of the asymptotic phase function 
evaluated at the limit cycle. For the iris system, the level curves appear to pinch 
together at the boundaries where they abruptly change direction, and this pinch- 
ing means the density of the level curves is highest at the square boundaries. One 
might naively expect that the phase response should be most sensitive to perturba- 
tions immediately adjacent to the saddle points; instead the peak sensitivity occurs 
at boundaries between regions dominated by one or another saddle. This efi^ect is 
qualitatively present both in the piecewise linear iris system and in the smooth sys- 
tem considered here. Within the piecewise linear regions of the iris system the flow is 
homogeneous; it is the boundaries that make the system nonlinear and it should not 
be surprising that that is where the greatest sensitivity occurs. In the sine system, 
it appears that near the saddles the flow is governed by the local linear approxima- 
tion, while in the region between saddles the flow is "more nonlinear" coinciding with 
greater sensitivity to perturbation. 

8.2. Comparison to the PRC near a homoclinic bifurcation. While the 
structure and asymptotic behavior of the phase response curve described in Theorem 
1 were derived specifically for the iris system, we expect the qualitative picture to 
hold for generic one-parameter families of limit cycles verging on a heteroclinic or 
a homoclinic bifurcation of vector fields as well. Treating the general case lies 
beyond the scope of this paper; however, it will be useful to compare our results with 
the analysis of phase resetting near a homoclinic bifurcation given in ([11], §3.1.3). 

In [11] the authors consider a vector field on M" with a hyperbolic saddle fixed 
point at the origin with a single unstable eigenvahie A„ and unstable eigenvalues A^j- 
with Au < Xs = minj[Asj]. For positive values of a bifurcation parameter /i, the 
system is assumed to have a limit cycle that spends the overwhelming majority of its 
time in a box B = [0, A]". Trajectories exit the box when the coordinate along the 
unstable eigendirection a; A and are instantaneously reinjected with x — e. When 
n = 2 the geometry is similar to that of the iris system (compare Figure 3 of [11] with 
our Figure 2.1b). The main difference is that in the homoclinic system, the value of 
the injection point along the unstable limit coordinate is assumed to be independent 
of the egress point, whereas in the iris system its dependence is specified by the map 
h = fe° fi (Figure 3.1). Equivalently, it is as if the derivative of the function h were 
zero at the fixed point v) rather than having finite slope. This situation would occur 
for the homoclinic system if the compression of the fiow taking trajectories from the 
egress boundary to the ingress boundary were sufficiently great. However, in analyzing 
the homoclinic system, one must also assume that the time spent along the portion of 
the limit cycle outside the box B is vanishingly small. Reconciling these two limiting 
processes for the general homoclinic case remains an interesting problem. 

For purposes of comparison, note that e/A in [11] corresponds to v) for the iris 
system, and the phase E [0, 2tt) in [11] corresponds to 2Tr(p in the iris system. Table 
8.1 compares the infinitesimal PRCs for the two systems. In [11] the perturbation 
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(u— As) log(l/u) 


{u—Xs) log(l/ti) 



Table 8.1: Phase response curve for homoclinic system [11] and the iris system, com- 
pared. For clarity, we write u, < u < 1, for = e/A, and the phase < <^ < 1 for 
both systems. Here s = u^; in the A ^ oo limit s ^ 0, and the systems coincide (up 
to a scaling factor Vx//S. related to the choice of perturbation direction). The iPRC 
for perturbations parallel to the unstable (resp. stable) eigenvector direction is given 
by Z„ (resp. Z^). 



is assumed to be in a particular direction corresponding to a voltage deflection in 
a neural model; the factor arises in the change of coordinates from the direction 
of voltage pcturbations relative to the unstable eigenvector direction. The flow is 
assumed to be infinitely compressive during the time between egress and reinjection, 
and consequently it is assumed that only components of the perturbation along the 
unstable cigendircction contribute to the iPRC. This assumption holds for the iris 
system in the limit as the unstable-to-stable eigenvalue ratio A„/As = A — >■ oo, but 
not for finite values of A. The A < oo case includes additional corrections; see Table 
8.1. 

Finally, assuming that the phase response is independent of displacement along 
the stable eigenvector direction is equivalent to assuming that the isochrons run par- 
allel to the stable eigenvector throughout the box B. For the iris system the slopes 
of the isochrons do approach zero as A — )■ oo, becoming more and more parallel to 
the stable eigenvector, but for A = oo the asymptotic phase and the isochrons are no 
longer well defined. For finite values of A the isochrons arc well defined, but are not 
parallel to the stable eigenvector. Compare the isochrons in the SW quadrant of the 
iris system shown in Figure 6.1. 

8.3. Phase resetting in the absence of an asymptotic phase. Determinis- 
tic limit cycles are not the only way to represent (approximately) periodic biological 
rhythms. Alternatives include heteroclinic or homoclinic attractors perturbed by 
small amplitude noise [8, 49] as well as fixed points of spiral sink type subject to 

small noisy perturbations [10, 31]. A system comprising a determinstic flow with a 
stable spiral fixed point, when perturbed by small to modest amounts of noise, will 
show noisy oscillatory trajectories that can be difficult to distinguish from a small 
amplitude limit cycle [10]. For example, the small oscillations in membrane potential 
of a nerve cell brought near to firing by a steady depolarizing current of modest size 
may be due either to "noisy spiral sink" dynamics or to "noisy small limit cycle" 
dynamics [46] . The classical definition of asymptotic phase breaks down for noisy sys- 
tems, for spiral sinks, and for heteroclinic or homoclinic orbits. Nevertheless "phase 
resetting" is actively studied in such systems both experimentally [21, 22, 46, 47] and 
theoretically [29, 40, 39]. For the Bonhocffer van dcr Pol equations, for example, 
Rabinovitcli and colleagues showed how to define an analog of the classical isochron, 
for point converging to a spiral sink along a distinguished trajectory (a "T-attractor") 
[40] . Generally speaking, noise efi'ects can expose the structure of bifurcations in neu- 
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ral systems through stochastic resonance [24] or spike time bifurcations [52, 53], and 
it would be useful to have a broader theory of "phase response" going beyond the 

case of deterministic limit cycles. 
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